El propósito de este artículo es visualizar los cambios en incidencia delictiva a raíz de la pandemia de COVID-19. Para ello, usaremos los datos de las carpetas de invertigacion de la PGJ de la CDMX, así como mapas a nivel AGEB y delegación.
Contenido:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import os
import geopandas as gpd
import descartes
import math
from datetime import datetime
from collections import Counter
from sklearn.cluster import DBSCAN
datos = pd.read_csv('inputs/carpetas-de-investigacion-pgj-de-la-ciudad-de-mexico.csv')
datos.head()
Para revisar la calidad de estos datos, revisaremos lo siguiente:
La unicidad por categoría: si los registros fueron hechos a mano se presta a que haya más de un registro por faltas de ortografía (por ejemplo, COYOACAN y COYOACÁN)
La cantidad de datos por categoría: se observa información hasta 1901, pero es probable que solo sea un hecho histórico registrado para ese año y no se tengan registros para años intermedios. Dado que tras observar la serie de tiempo se nota un incremento sustancial exactamente al iniciar 2016, lo cual se debe deber a un tema de registros, sólo consrvaremos la información disponible a partir de ese momento.
Convertimos las fechas a formato fecha. Para los valores faltantes tomamos la fecha de la variable fecha_inicio
Sólo conservamos datos de delitos ocurridos en la CDMX
# limpieza general
################################################################################
# solo conservamos registros de las top 16 alcaldias, que son las delegaciones, pues el resto son alcaldias menores del area metropolitana
ranking_delitos_alcaldia=datos[['alcaldia_hechos','delito']]\
.groupby('alcaldia_hechos').count()\
.reset_index()\
.sort_values(by='delito',ascending=False)
datos=datos[ datos['alcaldia_hechos'].isin( ranking_delitos_alcaldia.loc[:16,'alcaldia_hechos'] ) ]\
.reset_index().drop(columns='index')
del ranking_delitos_alcaldia
################################################################################
# limpiamos las fechas
datos['fecha']=datetime.strptime('2000-12-31 10:10','%Y-%m-%d %H:%M').date()
datos['aux']=datos['fecha_hechos'].apply(lambda x: str(x)[2])
aux_num=list(map(str,range(10)))
datos.loc[datos['aux'].isin(aux_num),'fecha']=datos.loc[datos['aux'].isin(aux_num),'fecha_hechos']\
.apply(lambda x: str(x)[:16])\
.str.replace('/','-').str.replace('T',' ')\
.apply(lambda x: datetime.strptime(x,'%Y-%m-%d %H:%M').date())
del aux_num
##########################
#n, se reemplazan por fecha_inicio que fluye sin problema
datos.loc[datos['aux']=='n','fecha']=datos.loc[datos['aux']=='n','fecha_inicio']\
.apply(lambda x: str(x)[:16])\
.str.replace('/','-').str.replace('T',' ')\
.apply(lambda x: datetime.strptime(x,'%Y-%m-%d %H:%M').date())
#/, solo tienen un formato de prsentación distinto
datos.loc[datos['aux']=='/','fecha']=datos.loc[datos['aux']=='/','fecha_hechos']\
.apply(lambda x: str(x)[:16])\
.str.replace('/','-').str.replace('T',' ')\
.apply(lambda x: datetime.strptime(x,'%d-%m-%Y %H:%M').date())
datos.loc[datos['aux']==' ','auxx']=datos.loc[datos['aux']==' ','fecha_inicio']\
.apply(lambda x: str(x)[2])
##########################
datos.loc[datos['auxx'].isin(['/','-']),'fecha']=datos.loc[datos['auxx'].isin(['/','-']),'fecha_inicio']\
.apply(lambda x: str(x)[:16]).str.replace('/','-').str.replace('T',' ')\
.apply(lambda x: datetime.strptime(x,'%d-%m-%Y %H:%M').date())
aux_num=list(map(str,range(10)))
datos.loc[datos['auxx'].isin(aux_num),'fecha']=datos.loc[datos['auxx'].isin(aux_num),'fecha_inicio']\
.apply(lambda x: str(x)[:16]).str.replace('/','-').str.replace('T',' ')\
.apply(lambda x: datetime.strptime(x,'%Y-%m-%d %H:%M').date())
del aux_num
datos=datos.drop(columns=['aux','auxx'])
################################################################################
# unicamente conservamos datos a partir de 2016
datos=datos[datos['fecha']>=datetime.strptime('01-01-2016','%d-%m-%Y').date()]
# generamos variables de mes_ano y semana
datos['mes_ano']=datos['fecha'].apply(lambda x: datetime.strptime((str(x.year))+'-'+str(x.month)+'-01','%Y-%m-%d').date())
datos['semana']=datos['fecha'].apply(lambda x: x.strftime("%U"))
Se propone comparar la incidencia de diferentes tipos de delitos antes y después de la llegada de la pandemia. Para ello, primero debemos determinar los conjuntos de información elegidos para cada período.
Para el período con COVID es bastante simple, pues sólo tenemos que ubicar el momento a partir del cuál la gente comenzó su confinamiento. Los datos muestran de manera clara la fecha en la que esto ocurrió: Abril. Casi todos los delitos disminuyen durante la primera mitad del año, debido a los efectos de la pandemia. Su posterior crecimiento y en algunos casos estabilización podrían considerarse como la "nueva normalidad".
datos.loc[ (datos['fecha']>datetime.strptime('01-01-2018','%d-%m-%Y').date()) &
(datos['fecha']<datetime.strptime('01-01-2021','%d-%m-%Y').date()),
['mes_ano','delito']]\
.groupby('mes_ano')\
.count()\
.plot(figsize=(15,4));
Para la determinación del período pre-COVID, lo más sencillo sería mantener todas las observaciones anteriores. Sin embargo, a pesar de que la información observada no muestra un choque releveante en la incidencia delictiva, han ocurrido cambios estructurales que obligan a hacer una exploración adicional.
En particular, valdría la pena buscar por cambios en patrones a raíz de la llegada de Claudia Sheinbaum a la jefatura de gobierno.
Desglosando los casos por tipo de delito podemos observar grandes cambios en sus respectivas categorías a lo largo del tiempo. Esto se puede deber a cambios en la estrategia de seguridad del gobierno, pero los cambios más grandes probablemente se deban a cambios en las tipificaciones de los delitos.
#datos.loc[datos['mes_ano']>=datetime.strptime('01-01-2020','%d-%m-%Y').date(),['mes_ano','categoria_delito','delito']]\
datos.loc[:,['mes_ano','categoria_delito','delito']]\
.groupby(['mes_ano', 'categoria_delito'])\
.count()\
.unstack()\
.plot(subplots=True,figsize=(10,60));
En la mayoría de los casos podemos observar cierta estabilidad, o al menos una tendencia sin grandes cambios, en el período a partir del cual llega Sheinbaum. Tomando en consideración posibles choques temporales derivados de la transición de gobierno, tomaremos la información 3 meses después de la toma de posesión de Sheinbaum, es decir, a aprtir de febrero de 2019.
datos=datos[datos['fecha']>=datetime.strptime('01-02-2019','%d-%m-%Y').date()].reset_index()
Con esta muestra reducida, vemos las correlaciones entre los distintos delitos, que sugieren algunas cosas:
delitos_por_colonia=datos[['categoria_delito','delito','colonia_hechos']]\
.groupby(['categoria_delito','colonia_hechos'])\
.count().reset_index()\
.pivot(index='colonia_hechos',columns='categoria_delito',values='delito')\
.fillna(0)
sns.heatmap(delitos_por_colonia.corr())
plt.show()
del delitos_por_colonia
Como análisis posterior se considera hacer pruebas de Perron de cambios de tendencia de series de tiempo.
Ahora, pasaremos a ver el cambio en incidencia delictiva per capita de cada tipo de crimen por delegación.
Se tomó la base de Indicadores de Peligro, Exposición y Vulnerabilidad a nivel AGEB del Atlas de Riesgo de la ciudad, con datos de población de 2010 (http://www.atlas.cdmx.gob.mx/datosabiertos.html). Posteriormente, se agregó la información al nivel de delegación, con una mapa de municipalidades de CONABIO (http://www.conabio.gob.mx/informacion/metadata/gis/muni_2018gw.xml?_httpcache%20=%20yes&_xsl=/db/metadata/xsl/fgdc_html.xsl&_indent%20=%20no) para el que conservamos únicamente el terriotorio correspondiente a la CDMX.
Para compararlos, tomaremos el mismo periodo de tiempo para cada año. Basándonos en los datos observados durante la pandemia, conservaremos la información de junio a noviembre, puesto que a partir de junio se "recuperó" la incidencia delictiva en la ciudad. De manera análoga, conservaremos los mismos meses para 2019 y los tomaremos como el periodo anterior.
#Por desgracia los datos de población no son muy estables, por lo que descartamos las AGEBs con menos de 600 habitantes para quitar observaciones sin información y con poca gente que tienen mucha varianza. Esto representa menos del 7% de las AGEBs, por lo que se conserva una parte sustancial del mapa.
datos['Periodo']='ninguno'
datos.loc[(datos['fecha']>=datetime.strptime('01-06-2019','%d-%m-%Y').date()) &
(datos['fecha']<datetime.strptime('01-12-2019','%d-%m-%Y').date()),'Periodo']='Pre'
datos.loc[(datos['fecha']>=datetime.strptime('01-06-2020','%d-%m-%Y').date()) &
(datos['fecha']<datetime.strptime('01-12-2020','%d-%m-%Y').date()),'Periodo']='Post'
datos=datos[datos['Periodo']!='ninguno']
datos_periodos=datos[['delito','categoria_delito','Periodo']].groupby(['categoria_delito','Periodo']).count().reset_index()
datos_periodos[datos_periodos['categoria_delito']!='DELITO DE BAJO IMPACTO']\
.pivot("categoria_delito", "Periodo", "delito")\
.plot(kind='barh',figsize=(12,7))
plt.show()
mapa = gpd.read_file("inputs/mapa_cdmx_ageb")
mapa['CVE_ALCALDIA']=mapa['CVEGEO'].str.slice(0,5)
mapa.head()
mapa_del = gpd.read_file("inputs/mapa_cdmx_del")
mapa_del = mapa_del[mapa_del['CVE_ENT']=='09']
mapa_del.head()
# genera base de puntos
points = gpd.GeoDataFrame(datos[['longitud','latitud','categoria_delito','delito','Periodo']],
geometry=gpd.points_from_xy(datos['longitud'],
datos['latitud']))
points=points[points.is_valid]
points.crs = mapa_del.crs
# total de delitos por tipo de delito y delegacion
total_delito_del=gpd.sjoin(points,mapa_del,
how='left',op='within')
total_delito_del=total_delito_del[['CVEGEO','categoria_delito','delito','Periodo']]\
.groupby(['CVEGEO','categoria_delito','Periodo'])\
.count().reset_index()
#poblacion por delegacion
pob_alcaldia=mapa[['CVE_ALCALDIA','ALCALDIA','POB_TOT']]\
.groupby(['CVE_ALCALDIA','ALCALDIA'])\
.sum().reset_index()
# al mapa de delegaciones sel añade la info total
mapa_del_plus=mapa_del\
.merge(total_delito_del,
how='left',on='CVEGEO')\
.merge(pob_alcaldia,
how='left',left_on='CVEGEO',right_on='CVE_ALCALDIA')
# tasa de crimen por zona
mapa_del_plus['tasa_delito']=mapa_del_plus['delito']/mapa_del_plus['POB_TOT']*100000
del points
del total_delito_del
del pob_alcaldia
mapa_del_plus.head()
datos_periodos_del=mapa_del_plus[['Periodo','categoria_delito','ALCALDIA','tasa_delito']]
datos_periodos_del.head()
intervalos_delitos=datos_periodos_del\
.groupby('categoria_delito')\
.agg({'tasa_delito':['min','max']})['tasa_delito'].reset_index()
intervalos_delitos.head()
A continuación se muestra una serie de mapas para observar la incidencia delicitiva por tipo de crimen y delegación. De este podemos sacar algunas conclusiones generales:
En particular, también nos dice algo de los efectos de la pandemia:
aux = mapa_del_plus\
.sort_values('Periodo',ascending=False)\
.groupby(['categoria_delito','Periodo'])
plt.figure(figsize=(10,64))
# itera por delitos y periodo
for i, ((aux_name,aux_per), aux_gdf) in enumerate(aux):
ax = plt.subplot(16, 2, i + 1)
aux_gdf.plot(ax=ax,column='tasa_delito',legend=True,
vmin=intervalos_delitos.loc[intervalos_delitos['categoria_delito']==aux_name,'min'],
vmax=intervalos_delitos.loc[intervalos_delitos['categoria_delito']==aux_name,'max'])
ax.set_title(aux_name+' - '+aux_per)
ax.set_aspect('equal', adjustable='datalim')
del aux
plt.show()
En los análisis subsecuentes, sólo trabajaremos con información de la "nueva normalidad"
datos_respaldo=datos.copy()
datos = datos[datos['Periodo']=='Post']
#datos_respaldo_2=datos.copy()
#datos=datos_respaldo.copy()
Por último, algo interesante sería ver, a nivel más desagregado, qué zonas son focos rojos actualmente. Una hipótesis podría ser que las zonas de aglomeración de actividades no escenciales que fueron cerradas por la pandemia cayeron en su índice de incidencia delictiva sustancialmente, mientras que los robos a casa habitación con violencia o similares aumentaron.
Para esto, mostraremos la ubicación de cada crimen en la ciudad y buscaremos construir conjuntos de los mismos mediante BDSCAN, una técnica de aprendizaje no supervisado.
points = gpd.GeoDataFrame(datos[['longitud','latitud','categoria_delito','delito','Periodo']],
geometry=gpd.points_from_xy(datos['longitud'],
datos['latitud']))
points=points[points.is_valid]
aux = points.groupby('categoria_delito')
plt.figure(figsize=(25,20))
for i, (aux_name, aux_gdf) in enumerate(aux):
ax = plt.subplot(4, 4, i + 1) # nrows, ncols, axes position
mapa_del.plot(ax=ax,color='white', edgecolor='grey')
aux_gdf.plot(ax=ax,markersize=3)
ax.set_title(aux_name)
ax.set_aspect('equal', adjustable='datalim')
del aux
plt.show()
datos[['categoria_delito','delito']]\
.groupby('categoria_delito').count().reset_index()
Como se puede observar, hay muy pocos registros de secuestros y robos en taxi, por lo que omitiremos dichas categorías para el siguiente ejercicio, así como los delitos de bajo impacto y los hechos no delictivos.
#categorias_delitos=list(set(datos['categoria_delito']))
categorias_delitos=['ROBO A PASAJERO A BORDO DEL METRO CON Y SIN VIOLENCIA',
'ROBO A TRANSEUNTE EN VÍA PÚBLICA CON Y SIN VIOLENCIA',
'ROBO A NEGOCIO CON VIOLENCIA',
'ROBO A CUENTAHABIENTE SALIENDO DEL CAJERO CON VIOLENCIA',
'LESIONES DOLOSAS POR DISPARO DE ARMA DE FUEGO',
'VIOLACIÓN',
'ROBO A REPARTIDOR CON Y SIN VIOLENCIA',
'ROBO A PASAJERO A BORDO DE MICROBUS CON Y SIN VIOLENCIA',
'ROBO DE VEHÍCULO CON Y SIN VIOLENCIA',
'HOMICIDIO DOLOSO',
'ROBO A CASA HABITACIÓN CON VIOLENCIA']
Para tener cierta perspectiva sobre las aglomeraciones, primero mostraremos un mapa de calor de población a nivel AGEB:
mapa.plot(column='POB_TOT',figsize=(11,9))
plt.show()
El algoritmo de BDSCAN toma como parámetros:
Para la distancia, tomamos 200 metros como la distancia adecuada para considerarlos dentro de una misma área. Para el número mínimo de observaciones, quizá con algunas objeciones válidas, se tomó el 2% del total de casos.
plt.figure(figsize=(25,80))
for i in range(len(categorias_delitos)):
# generación de puntos
points = gpd.GeoDataFrame(datos[(datos['Periodo']=='Post')&(datos['categoria_delito']==categorias_delitos[i])][['longitud','latitud','categoria_delito','delito']],
geometry=gpd.points_from_xy(datos.loc[(datos['Periodo']=='Post')&(datos['categoria_delito']==categorias_delitos[i]),'longitud'],
datos.loc[(datos['Periodo']=='Post')&(datos['categoria_delito']==categorias_delitos[i]),'latitud']))
points=points[points.is_valid]
# algoritmo de clasificación
clustering = DBSCAN(eps=.0115,min_samples=np.round(points.shape[0]*.02)).fit(np.array(points[['longitud','latitud']]))
ax = plt.subplot(6, 2, i + 1) # nrows, ncols, axes position
mapa.plot(ax=ax,color='white', edgecolor='grey')
points.plot(ax=ax,markersize=15,column=clustering.labels_,categorical=True)
ax.set_title(categorias_delitos[i])
ax.set_aspect('equal', adjustable='datalim')
#plt.tight_layout()
plt.show()
De esta visualización podemos sacar algunos comentarios: